clc
clear
close all
bar=waitbar(0,'建模中...');
%% Input initial parameters

InitialParameter = inputEssentialParameter();
InitialParameter = inputIntermediateBearing(InitialParameter );


% If you need change some parameters, please change the data in the struct:
% InitialParameter, then use establishModel( ) to get the different models

%% Establish models

% grid{1} = [1 4 1]; manualGrid{2} = [3]; 
grid = 'low';
Parameter = establishModel(InitialParameter,...
                           'gridfineness', grid,...
                           'isPlotModel',  true,...
                           'isPlotMesh',   true);
save('modelParameter','Parameter')  


%% Input calculate parameters
vStart = 10; % rad/s
vEnd = 1500; % rad/s
vNum = 1000;
accelerationTime = 3; % s
stableTime = 4; % s
transientPeriodNum = 100;
findTime = 2; % find the maximum value in the last 2 seconds
vmax = linspace(vStart, vEnd, vNum);
acceleration = vmax / accelerationTime;
transientTime = transientPeriodNum*2*pi ./ vmax; % s
if vStart == 0
    transientTime(1) = 0;
end

%% calcualte
% creat matrix save the result;
amplitude = zeros(Parameter.Mesh.dofNum, vNum); 

% waite bar
str = '开始计算...';
waitbar(0,bar,str)
len = vNum;

for iV = 1:1:vNum
    % change running statues
    Parameter.Status.vmax = vmax(iV);
    Parameter.Status.acceleration = acceleration(iV);

    % Generate the dynamic equation
    generateDynamicEquation(Parameter);    

    % calculate
    TSTART = 0;
    TEND = accelerationTime + transientTime(iV) + stableTime;
    SAMPLINGFREQUENCY = 10000;
    ISPLOTSTATUS = true;
    REDUCEINTERVAL = 1;
    tic
    [q, dq, t, convergenceStr] = calculateResponse(Parameter, [TSTART, TEND], SAMPLINGFREQUENCY,ISPLOTSTATUS,REDUCEINTERVAL);
    toc

    if ~isempty(convergenceStr)
        fprintf('%s \n', convergenceStr)
    end % end if 
    
    % find the maximum value
    amplitude(:, iV) = max(q(:, end-findTime*SAMPLINGFREQUENCY : end), [], 2); % return the maximum of each row 
    
    % waite bar
    str=['计算中...',num2str(100*iV/len),'%'];
    waitbar(iV/len,bar,str)
    
end % end for iV

% save the result
save('amplitudeFre', 'amplitude', 'vmax')

% close wait bar
delete(bar)

%% Plot

%name the label
dofNum          = Parameter.Mesh.dofNum;
nodeNum         = Parameter.Mesh.nodeNum;
dofOnNodeNo     = Parameter.Mesh.dofOnNodeNo;
figureIdentity  = cell(1,dofNum);
nodeNo          = 1;
dofInThisNode   = 0;
for iDof = 1:1:dofNum
    if nodeNo == dofOnNodeNo(iDof)
        dofInThisNode = dofInThisNode + 1;
    else
        nodeNo = nodeNo + 1;
        dofInThisNode = 1;
    end
    figureIdentity{iDof}=['Node-',num2str(dofOnNodeNo(iDof)),'-DOF-',num2str(dofInThisNode)];
end

for idof = 1:1:Parameter.Mesh.dofNum
    h = figure('visible', 'off', 'name', figureIdentity{idof});
    plot(vmax, amplitude(idof, :), 'linewidth',1.5,'color',[0.40784,0.5804,0.651]);
    set(gca, ...
        'Box'         , 'on'                        , ...
        'LooseInset'  , [0,0,0,0]                   , ...
        'TickDir'     , 'in'                        , ...
        'XMinorTick'  , 'off'                       , ...
        'YMinorTick'  , 'off'                       , ...
        'TickLength'  , [.01 .01]                   , ...
        'LineWidth'   , 0.5                         , ...
        'XGrid'       , 'on'                        , ...
        'YGrid'       , 'on'                        , ...
        'FontSize'    , 7                           , ... 
        'FontName'    , 'Times New Roman'           ,...
        'xtick'       , linspace(0,1500,11)) 
    xlabel('$\dot{\Phi}$ (rad/s)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
    isBearing = idof>= 65 ;
    isTranslation = rem(idof, 4)==1 || rem(idof, 4)==2;
    if isBearing
        ylabel('$q_{\mathrm{max}}$ (m)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
    elseif isTranslation
        ylabel('$q_{\mathrm{max}}$ (m)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
    else
        ylabel('$q_{\mathrm{max}}$ (rad)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
    end
    set(gcf,'Units','centimeters','Position',[6 6 7.2 4.5])
    set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
    figureName =  figureIdentity{idof};
    figurePath = 'G:/大学硕士/毕业论文/论文/result/fullModel/幅频特性/';
    savefig(h,[figurePath, figureName, '.fig']);
    print(h, [figurePath, figureName], '-depsc2')
    print(h, [figurePath, figureName], '-dpng')
end


%% mix some figure: low pressure rotor 
mycolor = [ 0.40784,0.5804,0.651;
            0.94902,0.37647,0.32157;
            191/255,155/255,111/255;
            242/255,203/255,5/255];
mixDofL = [1, 17, 33];
h = figure('visible', 'off', 'name', 'Low pressure rotor translation');
for idof = 1:1:length(mixDofL)
    plot(vmax, amplitude(mixDofL(idof), :), 'linewidth',0.5,'color',mycolor(idof,:)); hold on
end
set(gca, ...
    'Box'         , 'on'                        , ...
    'LooseInset'  , [0,0,0,0]                   , ...
    'TickDir'     , 'in'                        , ...
    'XMinorTick'  , 'off'                       , ...
    'YMinorTick'  , 'off'                       , ...
    'TickLength'  , [.01 .01]                   , ...
    'LineWidth'   , 0.5                         , ...
    'XGrid'       , 'on'                        , ...
    'YGrid'       , 'on'                        , ...
    'FontSize'    , 7                           , ... 
    'FontName'    , 'Times New Roman'           ,...
    'xtick'       , linspace(0,1500,11)) 
legend( '风扇','低压压气机','低压涡轮',...
        'Fontname', '宋体',...
        'FontSize',7,...
        'location', 'northwest');
xlabel('$\dot{\Phi}$ (rad/s)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
ylabel('$q_{\mathrm{max}}$ (m)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
set(gcf,'Units','centimeters','Position',[6 6 7.2 5.5])
set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
figureName =  'LPTranslation';
figurePath = 'G:/大学硕士/毕业论文/论文/result/fullModel/幅频特性/';
savefig(h,[figurePath, figureName, '.fig']);
print(h, [figurePath, figureName], '-depsc2')
print(h, [figurePath, figureName], '-dpng')

h = figure('visible', 'off', 'name', 'Low pressure rotor translation');
for idof = 1:1:length(mixDofL)
    plot(vmax, amplitude(mixDofL(idof)+2, :), 'linewidth',0.5,'color',mycolor(idof,:)); hold on
end
set(gca, ...
    'Box'         , 'on'                        , ...
    'LooseInset'  , [0,0,0,0]                   , ...
    'TickDir'     , 'in'                        , ...
    'XMinorTick'  , 'off'                       , ...
    'YMinorTick'  , 'off'                       , ...
    'TickLength'  , [.01 .01]                   , ...
    'LineWidth'   , 0.5                         , ...
    'XGrid'       , 'on'                        , ...
    'YGrid'       , 'on'                        , ...
    'FontSize'    , 7                           , ... 
    'FontName'    , 'Times New Roman'           ,...
    'xtick'       , linspace(0,1500,11)) 
legend( '风扇','低压压气机','低压涡轮',...
        'Fontname', '宋体',...
        'FontSize',7,...
        'location', 'northeast');
xlabel('$\dot{\Phi}$ (rad/s)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
ylabel('$q_{\mathrm{max}}$ (rad)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
set(gcf,'Units','centimeters','Position',[6 6 7.2 5.5])
set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
figureName =  'LPRotation';
figurePath = 'G:/大学硕士/毕业论文/论文/result/fullModel/幅频特性/';
savefig(h,[figurePath, figureName, '.fig']);
print(h, [figurePath, figureName], '-depsc2')
print(h, [figurePath, figureName], '-dpng')


%% mix some figure: high pressure rotor
mixDofH = [37, 53, 61];
h = figure('visible', 'off', 'name', 'High pressure rotor translation');
for idof = 1:1:length(mixDofH)
    plot(vmax, amplitude(mixDofH(idof), :), 'linewidth',0.5,'color',mycolor(idof,:)); hold on
end
set(gca, ...
    'Box'         , 'on'                        , ...
    'LooseInset'  , [0,0,0,0]                   , ...
    'TickDir'     , 'in'                        , ...
    'XMinorTick'  , 'off'                       , ...
    'YMinorTick'  , 'off'                       , ...
    'TickLength'  , [.01 .01]                   , ...
    'LineWidth'   , 0.5                         , ...
    'XGrid'       , 'on'                        , ...
    'YGrid'       , 'on'                        , ...
    'FontSize'    , 7                           , ... 
    'FontName'    , 'Times New Roman'           ,...
    'xtick'       , linspace(0,1500,11)) 
legend( '\fontname{宋体}高压压气机\fontname{Times New Roman}1',...
        '\fontname{宋体}高压压气机\fontname{Times New Roman}5',...
        '\fontname{宋体}高压涡轮',...
        'FontSize',7,...
        'location', 'northeast');
xlabel('$\dot{\Phi}$ (rad/s)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
ylabel('$q_{\mathrm{max}}$ (m)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
set(gcf,'Units','centimeters','Position',[6 6 7.2 5.5])
set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
figureName =  'HPTranslation';
figurePath = 'G:/大学硕士/毕业论文/论文/result/fullModel/幅频特性/';
savefig(h,[figurePath, figureName, '.fig']);
print(h, [figurePath, figureName], '-depsc2')
print(h, [figurePath, figureName], '-dpng')

h = figure('visible', 'off', 'name', 'High pressure rotor translation');
for idof = 1:1:length(mixDofH)
    plot(vmax, amplitude(mixDofH(idof)+2, :), 'linewidth',0.5,'color',mycolor(idof,:)); hold on
end
set(gca, ...
    'Box'         , 'on'                        , ...
    'LooseInset'  , [0,0,0,0]                   , ...
    'TickDir'     , 'in'                        , ...
    'XMinorTick'  , 'off'                       , ...
    'YMinorTick'  , 'off'                       , ...
    'TickLength'  , [.01 .01]                   , ...
    'LineWidth'   , 0.5                         , ...
    'XGrid'       , 'on'                        , ...
    'YGrid'       , 'on'                        , ...
    'FontSize'    , 7                           , ... 
    'FontName'    , 'Times New Roman'           ,...
    'xtick'       , linspace(0,1500,11)) 
legend( '\fontname{宋体}高压压气机\fontname{Times New Roman}1',...
        '\fontname{宋体}高压压气机\fontname{Times New Roman}5',...
        '\fontname{宋体}高压涡轮',...
        'FontSize',7,...
        'location', 'northeast');
xlabel('$\dot{\Phi}$ (rad/s)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
ylabel('$q_{\mathrm{max}}$ (rad)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
set(gcf,'Units','centimeters','Position',[6 6 7.2 5.5])
set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
figureName =  'HPRotation';
figurePath = 'G:/大学硕士/毕业论文/论文/result/fullModel/幅频特性/';
savefig(h,[figurePath, figureName, '.fig']);
print(h, [figurePath, figureName], '-depsc2')
print(h, [figurePath, figureName], '-dpng')


%% mix some figure: intermediate bearing node
mixDofI = [31, 61];
h = figure('visible', 'off', 'name', 'Intermediate bearing translation');
for idof = 1:1:length(mixDofI)
    plot(vmax, amplitude(mixDofI(idof), :), 'linewidth',0.5,'color',mycolor(idof,:)); hold on
end
set(gca, ...
    'Box'         , 'on'                        , ...
    'LooseInset'  , [0,0,0,0]                   , ...
    'TickDir'     , 'in'                        , ...
    'XMinorTick'  , 'off'                       , ...
    'YMinorTick'  , 'off'                       , ...
    'TickLength'  , [.01 .01]                   , ...
    'LineWidth'   , 0.5                         , ...
    'XGrid'       , 'on'                        , ...
    'YGrid'       , 'on'                        , ...
    'FontSize'    , 7                           , ... 
    'FontName'    , 'Times New Roman'           ,...
    'xtick'       , linspace(0,1500,11)) 
legend( '\fontname{宋体}节点\fontname{Times New Roman}6',...
        '\fontname{宋体}节点\fontname{Times New Roman}16',...
        'FontSize',7,...
        'location', 'northwest');
xlabel('$\dot{\Phi}$ (rad/s)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
ylabel('$q_{\mathrm{max}}$ (m)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
set(gcf,'Units','centimeters','Position',[6 6 7.2 5.5])
set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
figureName =  'IBTranslation';
figurePath = 'G:/大学硕士/毕业论文/论文/result/fullModel/幅频特性/';
savefig(h,[figurePath, figureName, '.fig']);
print(h, [figurePath, figureName], '-depsc2')
print(h, [figurePath, figureName], '-dpng')

h = figure('visible', 'off', 'name', 'Intermediate Bearing translation');
for idof = 1:1:length(mixDofI)
    plot(vmax, amplitude(mixDofI(idof)+2, :), 'linewidth',0.5,'color',mycolor(idof,:)); hold on
end
set(gca, ...
    'Box'         , 'on'                        , ...
    'LooseInset'  , [0,0,0,0]                   , ...
    'TickDir'     , 'in'                        , ...
    'XMinorTick'  , 'off'                       , ...
    'YMinorTick'  , 'off'                       , ...
    'TickLength'  , [.01 .01]                   , ...
    'LineWidth'   , 0.5                         , ...
    'XGrid'       , 'on'                        , ...
    'YGrid'       , 'on'                        , ...
    'FontSize'    , 7                           , ... 
    'FontName'    , 'Times New Roman'           ,...
    'xtick'       , linspace(0,1500,11)) 
legend( '\fontname{宋体}节点\fontname{Times New Roman}6',...
        '\fontname{宋体}节点\fontname{Times New Roman}16',...
        'FontSize',7,...
        'location', 'northwest');
xlabel('$\dot{\Phi}$ (rad/s)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
ylabel('$q_{\mathrm{max}}$ (rad)','Fontname', 'Times New Roman', 'FontSize',9,'interpreter','latex');
set(gcf,'Units','centimeters','Position',[6 6 7.2 5.5])
set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
figureName =  'IBRotation';
figurePath = 'G:/大学硕士/毕业论文/论文/result/fullModel/幅频特性/';
savefig(h,[figurePath, figureName, '.fig']);
print(h, [figurePath, figureName], '-depsc2')
print(h, [figurePath, figureName], '-dpng')










